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Dynamic  turbulence  modelling  in  large-eddy 
simulations  of  the  cloud-topped  atmospheric 

boundary  layer 

By  M.  P.  Kirkpatrick,  N.  N.  Mansour,  A.  S.  Ackerman  f,  and  D.  E.  Stevens  J 


1.  Motivation  and  Objectives 

The  use  of  large  eddy  simulation,  or  LES,  to  study  the  atmospheric  boundary  layer 
dates  back  to  the  early  1970s  when  DeardorflF  (1972)  used  a three-dimensional  simulation 
to  determine  velocity  and  temperature  scales  in  the  convective  boundary  layer.  In  1974 
he  applied  LES  to  the  problem  of  mixing  layer  entrainment  (Deardorff  1974)  and  in  1980 
to  the  cloud-topped  boundary  layer  (Deardorff  19806).  Since  that  time  the  LES  approach 
has  been  applied  to  atmospheric  boundary  layer  problems  by  numerous  authors. 

While  LES  has  been  shown  to  be  relatively  robust  for  simple  cases  such  as  a clear, 
convective  boundary  layer  (Mason  1989),  simulation  of  the  cloud-topped  boundary  layer 
has  proved  more  of  a challenge.  The  combination  of  small  length  scales  and  anisotropic 
turbulence  coupled  with  cloud  microphysics  and  radiation  effects  places  a heavy  bur- 
den on  the  turbulence  model,  especially  in  the  cloud-top  region.  Consequently,  over  the 
past  few  decades  considerable  effort  has  been  devoted  to  developing  turbulence  models 
that  are  better  able  to  parameterise  these  processes.  Much  of  this  work  has  involved 
taking  parameterisations  developed  for  neutral  boundary  layers  and  deriving  corrections 
to  account  for  buoyancy  effects  associated  with  the  background  stratification  and  local 
buoyancy  sources  due  to  radiative  and  latent  heat  transfer  within  the  cloud  (see  Lilly 
1962;  Deardorff  1980o;  Mason  1989;  MacVean  & Mason  1990,  for  example).  In  this  pa- 
per we  hope  to  contribute  to  this  effort  by  presenting  a number  of  turbulence  models 
in  which  the  model  coefficients  are  calculated  dynamically  during  the  simulation  rather 
than  being  prescribed  a priori. 

The  dynamic  procedure,  first  proposed  by  Germano  et  al.  (1991),  computes  model 
coefficients  using  information  contained  in  the  resolved  velocity  and  scalar  fields.  In  this 
sense  dynamic  models  can  be  considered  self-caUbrating,  a feature  that  makes  them  an 
appealing  choice  for  dealing  with  the  complex  interactions  between  the  hydrodynamics, 
radiation  and  cloud  microphysics  occurring  within  clouds.  Nevertheless,  while  dynamic 
turbulence  models  have  been  used  with  considerable  success  for  complex  engineering 
flows  (see  Boivin  et  al.  2000;  Branley  & Jones  2001,  for  example),  their  application  to 
atmospheric  flows  has  been  very  limited. 

The  objective  of  this  paper  is  to  present  a set  of  dynamic  turbulence  models  written 
in  a form  appropriate  for  large  eddy  simulations  of  the  atmospheric  boundary  layer  in 
which  the  flow  is  described  using  equations  of  motion  in  anelastic  form  . The  models  are 
tested  using  simulations  of  a nocturnal  marine  stratocumulus  cloud  observed  during  dur- 
ing the  second  D3mamics  and  Chemistry  of  Marine  Stratocumulus  (DYCOMS-II)  field 
experiment.  This  test  case  includes  a number  of  features  that  typically  cause  problems 
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in  simulations  of  the  atmospheric  boundary  layer,  including  the  presence  of  a strong 
temperature  inversion  at  the  level  of  the  cloud  top,  and  positive  feedback  loops  involv- 
ing turbulent  entrainment  across  the  inversion,  radiation  and  cloud  microphysics.  As 
such,  these  simulations  provide  a good  test  for  the  utility  of  our  models.  Simulations 
are  performed  using  the  atmospheric  boundary  layer  LES  code,  DHARMA,  (Stevens  & 
Bretherton  1999;  Stevens  et  al.  2000),  with  dynamic  turbulence  models  computed  using 
the  new  LLAMA  code  presented  here. 


2.  Dynamic  turbulence  models  for  the  anelastic  equations 

The  dynamics  of  the  cloud-topped  atmospheric  boundary  layer  can  be  described  using 
equations  for  conservation  of  mass,  momentum,  liquid  water  potential  temperature  and 
total  water  mixing  ratio.  These  equations  are  written  in  the  anelastic  form  of  Ogura  & 
Phillips  (1962)  in  which  the  thermodynamic  variables  such  as  pressure  p are  decomposed 
into  an  isentropic  base  state  po  (corresponding  to  a uniform  potential  temperature  6o) 
and  a dynamic  component.  Following  Clark  (1979),  the  dynamic  component  is  further 
decomposed  into  an  initial  environmental  deviation  in  hydrostatic  balance  pi  and  a time- 
evolving  d5mamic  perturbation  p2  to  give 


p{x,  y,  2,  t)  = po{z)  + Pi  (z)  + p2(x,  y,  z,  t). 


(2.1) 


After  subtracting  hydrostatic  balances,  the  resulting  continuous  equations  are  written. 


dpoUj  digpUjUj) 

dt  dxj 

dQoet  d(Qoetuj) 

dt  dxj 

dgoqt  djgoqtUj) 

dt  dxj 

djgoUj) 

dXj 


dP2  , r ^Q0^v2  „ , , TT 

— dx-  ' 0 ^ijkQofj'^k  "K 

= He;, 

= 0. 


(2.2) 

(2.3) 

(2.4) 

(2.5) 


Here  Uj  are  the  Cartesian  components  of  the  velocity  vector,  g the  density,  Sij  the  Kro- 
necker  delta,  eyfe  the  permutation  tensor,  fj  the  Coriolis  parameter,  g the  acceleration 
due  to  gravity,  qt  the  total  water  mixing  ratio  and  9[  = {6i  — Bp) /Op  a scaled  liquid  wa- 
ter potential  temperature.  Total  water  mixing  ratio  is  the  sum  of  the  liquid  and  vapour 
mixing  ratios, 

qt=qc  + Qv  = , (2.6) 

6d 

where  gc,  g^  and  gd  are  the  density  of  the  condensed  water,  water  vapour  and  dry  air 
respectively.  Liquid  water  potential  temperature  is  defined  as 


6i=e- 


Lqc 


(2.7) 


Here  L is  the  latent  heat  of  vaporisation,  Cpd  the  specific  heat  at  constant  pressure  for 

dry  air  and  ttq  = ( ^ with  Pref  a reference  pressure  and  Rd  the  gas  constant  of 

dry  air.  The  virtual  potential  temperature  appearing  in  the  buoyancy  term  of  the 
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momentum  equations  is  given  by 


Sv  — 6 do 


(2.8) 


where  Ry  is  the  gas  constant  for  water  vapour.  The  terms  Hy.,  Hg-  and  are  source 
terms  that  include  parameterisations  for  physical  processes  such  as  radiation,  precipita- 
tion and  the  effects  of  the  global  atmospheric  circulation. 

In  a large  eddy  simulation,  the  equations  are  filtered  to  remove  from  the  solution  those 
fluctuations  that  cannot  be  resolved  by  the  numerical  method.  For  the  anelastic  equations 
it  is  convenient  to  use  a density-weighted  or  Favre  filter,  where  a Favre  filtered  variable 
is  defined  as  (f>  = q4>/q-  Application  of  this  filter  to  the  equations  gives 


dpoUi  , d{QoUiUj)  dp2  , r _Qo^v2  _ ^ , TT 

~dT  + arr,.  ■ ^ 


d{QoO*iUj) 

dt  dxj 

dpoQt  djQoqtUj) 
dt  dXj 

9(.QoUj) 
dec  A 


Hbt- 


dxj  ’ 
dxj  ’ 

0, 


with  subgrid  scale  stresses  and  fluxes  given  by 

'^ij  ” 00  » 

) 


ye;  = 00 

Iqt  = ?o  - QtUj) . 


9Tjj 
dxj  ’ 


(2.9) 

(2.10) 

(2.11) 

(2.12) 


(2.13) 

(2.14) 

(2.15) 


We  have  assumed  here  that  fj  is  constant  so  that  the  Coriolis  term  is  linear  - an  assump- 
tion that  is  valid  when  the  domain  is  small  compared  to  the  Rossby  radius  of  deformation. 
The  source  terms  H are  shown  without  an  overbar  since  they  are  parameterisations  rather 
than  exact  terms. 

The  Smagorinsky  model  (Smagorinsky  1963)  for  the  subgrid  scale  stress  is 


'^ij  — ^^ij'^kk  — ‘^do^mDijCsi 

where  Dij  is  the  strain  rate  tensor, 


(2.16) 




2 \dxj  dxi)  3 dxk' 


(2.17) 


Here  and  elsewhere  in  this  paper  the  superscript  “ is  used  to  denote  the  anisotropic  part 
of  a tensor.  The  eddy  viscosity  Km  is  given  by 


Km  = CA^ 


ID 


(2.18) 


with 


D\ 


2DijDij  and  C a dimensionless  coefficient.  Stratification  effects  are  ac- 
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counted  for  using  the  correction  of  Lilly  (1962), 


Cb  = (1 

Cb  = 0 


Pisgs  ^ 1 


where  the  subgrid  scale  flux  Richardson  number  is  estimated  using 


P^sgs  — 


g/00  dO^jdz 


(2.19) 

(2.20) 


(2.21) 


Prags  is  the  subgrid  scale  turbulent  Prandtl  number,  which  we  set  to  a constant  value  of 
Prags  = 0.4. 

The  mixed  model  of  Bardina  et  al.  (1983),  combines  the  Smagorinsky  model  with  the 
scale  similarity  model  of  Bardina  et  al.  (1980).  The  model  is  derived  by  rewriting  the 
subgrid  scale  stress  tensor  in  terms  of  resolved  and  unresolved  parts  of  the  velocity  field 
m =Ui+u[  using  the  decomposition  of  Germane  (1986).  For  the  anelastic  equations  this 
decomposition  takes  the  form 

7v^  = LL  + C^  + R^,  (2.22) 

where  the  modified  Leonard,  cross  and  SGS  Re30iolds  terms  are 

L^j  — pQ^UiUj  UiUj)  (2.23) 

C*j  = Qoiuiu'j  + u'fi.  - UiU'j  - u'iUj)  (2.24) 

P*j  = Qoiu'iU'^  - u'iu'j).  (2.25) 

The  modified  Leonard  term,  which  contains  only  resolved  scale  quantities,  can  be  calcu- 
lated explicitly  leaving  the  modified  cross  and  SGS  Reynolds  terms  to  be  parameterised 
using  the  Smagorinsky  model. 


’’’ij  ~~  ^ij 


2qoCA'^  D DijCB. 


(2.26) 


Martin  et  al.  (2000)  note  that,  with  the  mixed  model,  the  isotropic  part  of  the  subgrid 
scale  stress  ry  can  be  represented  at  least  partially  by  the  isotropic  part  of  the  modified 
Leonard  term.  In  fact  they  obtained  the  closest  agreement  with  experiment  for  their 
compressible  flow  test  cases  when  they  assumed  that  the  isotropic  part  of  L*j  models  the 
whole  of  \SijTkk-  We  have  adopted  the  same  approach  in  Eq.(2.26). 

In  dynamic  versions  of  these  models,  the  model  coefficient  C is  calculated  dynamically. 
To  achieve  this,  a test  filter  is  applied  to  the  velocity  and  scalar  fields  to  extract  informal 
tion  from  the  smallest  resolved  scales.  Application  of  a spatial  test  filter  (denoted  here 
by  a caret)  to  the  filtered  momentum  equations  gives 


&QoUi  d(goUiUj) 

r\  • r\ 


^P2  , s „6o^v2 


^ijkQo'^kfj  “b  Hui 


(2.27) 


In  order  to  rewrite  this  equation  in  a form  similar  to  Eq.(2.9)  we  adopt  a Favre  test  filter 
^ = 94>/q  giving 


dQoUiUj)  dpl  , j 'qo0v2  ^ 

dxg  ~ Oo  + Hu,  . 


(2.28) 


Since  density  varies  only  in  the  z direction,  this  equation  can  be  simplified  considerably 


LES  of  the  cloud-topped  atmospheric  boundary  layer 


27 


by  using  two-dimensional  horizontal  test  filter  rather  than  a three-dimensional  one.  The 
Favre  test  filter  is  then  equivalent  to  a spatial  test  filter, 


= 

and  the  density  field  is  unchanged  by  the  test  filtering  operation,  that  is, 

^0  ~ Qo- 

Prom  Eqs.(2.29)  and  (2.30),  with  horizontal  test  filtering  Eq.(2.28)  becomes 


(2.29) 

(2.30) 

(2.31) 

(2.32) 

Alternatively,  application  of  the  grid  and  test  filters  together  to  the  continuous  equations 
gives 


dpi  ^ ~ , TT  9?ij 

- + ~ ^ijkQoUkfj  + Hu,  - 


dt 


dxi 


Extracting  the  resolved  stress  = gQ(uiUj  — UiUj)  then  gives 


dgpUi  d{gpUiUj) 

dt  dxj 


dp*  ^ 

— = + Sisg^g  - ^ijkQo^kfj  + Hu, 


dTjj  _ diiij 
dxi  dxi 


dgpUi  , d(goUiUj)  dpX  ^ , ^ , rr 

dt  dXj  ~ (ijkQoUkfj  + Hu, 


where 


(2.33) 

(2.34) 

(2.35) 

The  dynamic  procedure  assumes  that  the  same  parameterisation  that  is  used  to  pa- 
rameterise  the  SGS  stress  at  the  grid  filter  level  can  be  used  to  model  the  test  level 
stress  Tij.  The  test  level  stress  for  the  Smagorinsky  and  mixed  models  are  written  as 


'^ij  — Qo  • 

Comparison  of  Eqs.(2.32)  and  (2.33)  gives  the  Germano  identity, 


£jij  — Tij  Tij . 


Tfj  = -2gpCA^ 


D 


DijCs, 


and 


D 


DijCs, 


Cb  — (l-—  Risgsj 


X 1/2 


Hisgs  ^ 1 


Tij  = L:J  -2gpCA\ 
respectively.  Here  Cq  is  given  by 

^1  — Ris 

Cb  = 0 Risgs  > 1) 

with  the  test  level  Richardson  number  given  by 

g/00  dOuldz 

■^"^sgs  — • 

Prsgs\D\^ 

The  test  level  modified  Leonard  term  in  the  mixed  model  LlJ  is  found  by  decomposing  Tj 


(2.36) 

(2.37) 

(2.38) 

(2.39) 

(2.40) 
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in  a manner  similar  to  the  decomposition  of  (see  Eq.(2.22)).  Substituting  Ui  = Ui+u" 
into  Eq.(2.34)  gives 

+ + (2.41) 

where  the  modified  Leonard,  cross  and  Reynolds  terms  at  the  test  level  are 


L*f  — Qo{uiUj  - UiUj)  (2.42) 

0*^  = goluiUj  + u"uj  - Uiu"  - u'luj)  (2.43) 

(2.44) 


In  Eqs.(2.36-2.44)  we  have  again  used  the  fact  that  we  assume  a horizontal  test  filter  so 
that  Eqs. (2.29-2.30)  can  be  used  to  simplify  the  equations. 

The  decomposition  shown  in  Eqs. (2.42-2.44)  is  different  from  that  derived  by  Zang 
et  al.  (1993),  who  used  a decomposition  based  on  Ui  =Ui+u'^  at  both  the  grid  and  test 
level,  rather  than  the  form  Ui  =Ui+  u"  used  here  for  the  test  level.  As  noted  by  Vreman 
et  al.  (1994),  the  latter  is  more  consistent  with  the  dynamic  procedure  since  the  test  level 
stress  should  be  written  in  terms  of  test  filtered  velocity  components  in  order  to  make 
the  test  level  parameterisations  analogous  to  those  used  at  the  grid  level. 

Substituting  the  parameterisations  for  the  subgrid  scale  stress  Tjj  (Eq.(2.26))  and  test 
level  stress  Tij  (Eq.(2.37))  into  the  Germano  identity  (Eq.(2.35))  gives 


QoiuiUj  - UiUj)  = eoi^iUj  - UiUj)  - 2goCA^ 


D 


DijCs 


(2.45) 


-Q^luiUj -UiUj) -\-2qqCA'^  D DijCB- 


Removing  the  constant  factor  Qq,  and  contracting  tensors  using  the  least  squares  approach 
of  Lilly  (1992)  gives  an  equation  for  the  model  coefficient  in  the  dynamic  mixed  model, 


= -- 


< Mij{L*j  Hij)  > 


< 2MkiMki  > 


where 


Aij  — UiUj  UiUjj 


Mij  = 


D 


DijCs  — 


D 


DijCs, 


(2.46) 

(2.47) 

(2.48) 


Hij  — {UiUj  — UiUj)  - {UiUj  — UiUj), 
a is  the  ratio  of  the  test  and  grid  filters 

a = A/A, 


(2.49) 

(2.50) 


and  < > indicates  averaging  on  horizontal  planes.  Equivalently,  for  the  dynamic  Smagorin- 
sky  model. 


< > 


CA^  = / . (2.51) 

< 2MkiMki  > ^ ' 

The  models  outlined  above  use  the  original  d3mamic  procedure  of  Germano  et  al. 
(1991)  and  require  averaging  in  homogeneous  directions  to  remain  stable.  Ghosal  et  al. 
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(1995)  recast  the  dynamic  procedure  as  a constrained  variational  problem  and  used  this 
approach  to  develop  a “dynamic  localization  model”  that  does  not  require  averaging. 
This  model,  however,  requires  the  solution  of  an  integral  equation  and  is  computationally 
expensive.  A less  expensive  alternative  proposed  by  Piomelli  & Liu  (1995)  involves  finding 
an  approximate  solution  to  the  integral  equation  by  using  the  value  of  C at  the  previous 
time  step  to  give  a first  approximation  C*.  Applying  their  procedure  gives  localised 
versions  of  the  models.  The  coeflBcient  for  the  localised  dynamic  Smagorinsky  model  is 
given  by 


C = - 


^■AkiAki 


(2.52) 


where. 


Aij  = 


Bij  = \d 


D 


DijCs, 


I DijCs- 

The  coefficient  for  the  localised  dynamic  mixed  model  is 

(£t.  - Hjj  - 2C^ij)Aij 


C = 


1 Aij  Aij 


(2.53) 

(2.54) 


(2.55) 


with  Aij  and  Bij  given  by  Eqs.(2.53)  and  (2.54)  respectively  and  Hij  given  by  Eq.(2.49). 

We  note  here  that  by  using  a horizontal  test  filter  we  have  reduced  the  complexity  of 
our  models  considerably.  With  three-dimensional  test  filtering  the  models  would  be  of 
similar  complexity  to  the  equivalent  dynamic  models  for  compressible  flow  whereas,  when 
a horizontal  filter  is  used,  the  models  for  the  anelastic  equations  become  similar  to  those 
for  incompressible  flow.  In  fact,  if  the  discretisation  operator  is  taken  as  being  equivalent 
to  the  first  level  of  grid  scale  Favre  filtering  (that  is,  we  assume  implicit  filtering  at  the 
grid  scale),  then  the  only  difference  in  the  implementation  occurs  in  the  mixed  model, 
where  the  second  grid  scale  filter  in  the  formula  for  Hij  (Eq.(2.49))  is  a Favre  filter  that 
must  be  applied  explicitly. 

An  interesting  feature  of  the  localised  dynamic  model  is  that,  if  unconstrained,  it  can 
calculate  negative  values  for  the  model  coefficient,  and  thus  predict  regions  of  negative 
eddy  viscosity.  While  such  regions  may  be  interpreted  as  regions  of  backscatter,  or  energy 
transfer  from  small  to  large  scales,  Ghosal  et  al.  (1995)  and  Carati  et  al.  (1995)  note  that 
the  dynamic  Smagorinsky  model  has  no  information  about  how  much  energy  is  actually 
available  in  the  subgrid  scales.  Consequently,  there  is  no  mechanism  by  which  such  a 
“counter-gradient  diffusion”  model  for  backscatter  can  be  turned  off  once  all  the  subgrid 
scale  energy  has  been  removed.  For  this  reason  we  chose  to  clip  the  eddy  viscosity  at 
zero.  For  the  localised  model  backscatter  is  therefore  included.  For  the  plane  averaged 
dynamic  Smagorinsky  model,  however,  the  effect  of  backscatter  is  included  in  an  average 
sense  through  a reduction  of  the  plane  averaged  model  coefficient.  As  noted  above,  in  the 
mixed  model  the  Leonard  term  allows  backscatter,  so  there  is  no  need  to  resort  to  using 
a negative  eddy  viscosity. 

The  Smagorinsky  and  mixed  models  for  the  subgrid  scale  fluxes  of  the  scalar  variables 
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d\  ^Cb, 
I dxj 


l<t>  = - 4>Uj)  - QqC^A^ 


D 


dXj 


(2.56) 

(2.57) 
a term  repre- 


\ I WJkj 

Here  is  a generic  scalar  representing  61  or  qt-  The  mixed  model  includes 
senting  the  scalar  equivalent  of  the  modified  Leonard  term, 

Equations  for  the  model  coefficients  are  then  written  in  a form  analogous  to  that  used 
for  the  subgrid  scale  stress.  For  the  plane  averaged  dynamic  Smagorinsky  model, 

^ p.p. 


(2.58) 


^ A 2 _ < FjEj  > 

where  the  resolved  flux  Ej  is  given  by 

Ej  = -$Uj, 

and 

Fj  = o? 

For  the  dynamic  mixed  model. 


D 


\D 


d% 

dXj 


Cb. 


- — < aa  > ' 


where 


Gj  = {(j)Uj  {^Uj  -^Uj). 

The  coefficient  for  the  localised  dynamic  Smagorinsky  model  is  gi 


given  by 


^ (E,-CiQl)Pl 

°* KK 


where 


Pj  = 

Qi 


D 


d<p 

dxj 


Cb. 


The  coefficient  for  the  localised  dynamic  mixed  model  i 


IS 


^ {Ej-Gi-ClQi)P^ 

KK 


(2.59) 

(2.60) 
(2.61) 
(2.62) 

(2.63) 

(2.64) 

(2.65) 

(2.66) 

(2.67) 


PkPk 

with  Pj  and  Qj  given  by  Eqs.(2.65)  and  (2.66)  respectively,  and  Gj  given  by  Eq.(2.63). 

Close  to  the  rough  surface  at  the  bottom  of  the  domain  the  size  of  the  turbulent  eddies 
decreases  and  the  turbulence  models  described  above  do  not  give  the  correct  dissipation. 
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The  standard  Smagorinsky  model  is  over-dissipative  while  the  dynamic  models  are  under- 
dissipative.  With  the  standard  Smagorinsky  model  we  use  a damping  function  suggested 
by  Mason  (1994), 


1 


1 

— + 


1 


A2  ' {k{z  + zo)y' 

where  k is  the  Von  Karman  constant  and  zq  the  roughness  length.  With  all  models  we 
also  add  a surface  layer  stress  term  suggested  by  Brown  et  al.  (2001).  This  stress  takes 
the  form 


Tgfc  — J ^sfc^iz)po\'^\4'^^f 


(2.69) 


where  Csfc  is  a scaling  factor,  a(z)  a smoothing  function,  |u|  the  horizontal  wind  speed 
and  (j)  the  dependent  variable.  At  present  we  follow  Chow  & Street  (2002)  and  choose 
Csfc  to  give  Tsfc  equal  to  half  the  surface  flux  at  the  first  grid  point  above  the  wall. 


3.  Simulation  results 

The  test  case  chosen  for  the  simulations  is  the  first  research  flight  (RFOl)  of  the  second 
Dynamics  and  Chemistry  of  Marine  Stratocumulus  (DYCOMS-II)  field  experiment.  The 
computational  domain  extends  3.2km  in  the  two  horizontal  directions  and  1.5km  in  the 
vertical  direction,  with  96  x 96  x 128  cells  in  the  x,  y and  z directions  respectively.  Cell- 
spacing  is  uniform  in  the  horizontal  directions  and  stretched  in  the  vertical  direction  to 
give  cells  of  height  5m  close  to  the  bottom  surface  and  in  the  vicinity  of  the  inversion. 
This  grid  is  typical  of  that  currently  used  in  large  eddy  simulations  of  the  atmospheric 
boundary  layer.  The  spatial  discretisation  uses  a second-order  centered  scheme  with 
no  flux  limiting,  while  time  integration  uses  a second  order  Runge-Kutta  scheme.  The 
simulation  is  integrated  over  a period  of  four  hours.  Further  details  of  the  simulation 
set-up  are  described  in  Kirkpatrick  et  al.  (2003)  and  will  not  be  repeated  here. 

In  the  following  we  compare  results  obtained  with  four  of  the  SGS  turbulence  models 
described  above:  the  standard  Smagorinsky  (SM),  dynamic  Smagorinsky  (DSM),  dy- 
namic mixed  (DMM)  and  localised  dynamic  Smagorinsky  model  (LDSM).  The  standard 
Smagorinsky  model  uses  a constant  coefficient  of  Cs  = C^/2  = q.18  and  a Prandtl 
number  of  Prggs  = 0.4.  For  reference,  we  also  show  results  obtained  with  no  SGS  model. 

Figure  1 shows  vertical  profiles  of  condensed  water,  liquid  water  potential  temperature, 
eddy  viscosity  and  diffusivities  and  SGS  Prandtl  number  predicted  by  the  various  models. 
These  profiles  are  calculated  by  averaging  over  horizontal  planes  as  well  as  over  a time 
period  of  30  minutes.  All  profiles  are  calculated  for  the  last  half  hour  of  the  simulation, 
that  is  for  t = 3.5  — 4h.  The  d3mamic  Smagorinsky  models  give  quantitatively  similar 
profiles  to  the  standard  Smagorinsky  model.  LDSM  gives  somewhat  higher  average  eddy 
viscosity  and  diffusivities  than  DSM.  This  is  because  negative  values  of  eddy  viscosity 
predicted  by  LDSM  are  clipped,  whereas  with  DSM  conditions  that  would  produce  lo- 
cal negative  values  have  the  effect  of  reducing  the  plane  average.  The  undipped  model 
coefficient  profiles  for  DSM  and  LDSM  (not  shown)  were  found  to  be  almost  identical. 
The  djmamic  mixed  model  gives  lower  values  of  eddy  viscosity  and  diffusivity  because 
the  Leonard  term  component  of  the  SGS  stress  is  included  separately.  All  models  predict 
a SGS  turbulent  Prandtl  number  close  to  the  standard  value  of  0.4  through  most  the 
boundary  layer,  although  DSM  predicts  an  increase  to  a value  of  approximately  2 in  the 
region  just  below  the  inversion.  The  reason  for  this  behaviour,  which  is  not  seen  with  the 
other  models,  is  not  clear  at  present. 
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Time  (h) 

Figure  2.  Time  series  of  inversion  height,  cloud  fraction  and  liquid  water  path: no  SGS 

model,  ■ ■ ■ SM, DSM, LDSM, DMM. 


Figure  2 shows  time  series  of  inversion  height  Zinv,  liquid  water  path  LWP  and  cloud 
fraction.  The  entrainment  rate  E can  be  estimated  from  the  rate  of  change  of  Zinv  using 
the  formula, 

_ dZinv  /n  t\ 

E — — Wsxtb-  (^'1) 

Here  Wsub  is  the  subsidence  velocity,  which  for  the  present  case  is  estimated  to  be 
— 0.3cms"^.  dzinv/dt  is  calculated  as  the  average  rate  of  change  over  the  period  t = 2— 4h. 
The  dynamic  Smagorinsky  models  give  an  entrainment  rate  of  approximately  0.36- 
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Figure  3.  Eddy  diffusivity,  buoyancy  frequency,  and  SGS  and  resolved  buoyancy  production 
at  cloud-top:  • • • SM, DSM,  — • — LDSM,  — • • • — DMM. 

0.37cms"^ , while  the  dynamic  mixed  model  gives  0.40cms~  ^ . These  results  are  in  excellent 
agreement  with  the  entrainment  rate  E = 0.38±0.04cms“^  estimated  from  observational 
measurements.  The  standard  Smagorinsky  model,  however,  overpredicts  entrainment  sig- 
nificantly giving  E w 0.55cms~^.  The  entrainment  rate  with  no  SGS  model  is  0.34cms“^, 
which  is  slightly  lower  than  that  obtained  with  the  dynamic  models,  but  still  within  the 
error  bounds  of  the  observational  measurements.  The  observed  liquid  water  path  and 
cloud  fraction  remained  approximately  constant  during  the  experiment.  This  is  also  well 
predicted  with  dynamic  models  and  no  SGS  model,  while  with  the  standard  Smagorinsky 
model  the  results  for  these  parameters  again  deviate  significantly. 

To  help  explain  these  differences.  Figure  3 shows  a comparison  of  the  eddy  diffusivity 
Kh,  squared  buoyancy  frequency  N'^,  resolved  buoyancy  production  B and  subgrid  scale 
buoyancy  production  Bsgs-  Subgrid  scale  buoyancy  production  is  approximated  as 

« -KhN\  (3.2) 

The  buoyancy  frequency  is  calculated  using  a method  described  by  Stevens  et  al.  (2002). 
This  approach  is  an  extension  of  the  work  of  MacVean  & Mason  (1990)  and  takes  into 
account  changes  in  stability  that  occur  when  subgrid  scale  mixing  causes  condensation  or 
evaporation.  Figure  3 shows  the  eddy  diffusivity  predicted  by  the  standard  Smagorinsky 
model  decreasing  through  the  inversion  as  a result  of  the  stable  stratification  which  drives 
the  buoyancy  correction  coefficient  Cg  in  the  model  to  zero  (see  Eq.(2.16)).  There  is, 
however,  a significant  overlap  region  in  which  both  the  eddy  diffusivity  and  buoyancy 
frequency  are  large.  Prom  Figure  3 and  Eq.(2.16)  it  is  clear  that  this  overlap  region 
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Figure  4.  Second  and  third  moments  of  vertical  velocity.  Simulation  results  (resolved  component 

only): no  SGS  model,  • • • SM, DSM, LDSM, DMM.  Measurements; 

radar  data  (open  circles),  in  situ  measurements  (filled  circles). 


also  constitutes  a zone  of  negative  subgrid  scale  buoyancy  production.  Meanwhile  the 
resolved  buoyancy  production  with  SM  is  similar  to  that  seen  with  the  other  models. 
Negative  buoyancy  production  is  associated  with  conversion  of  mechanical  energy  into 
potential  energy.  In  the  present  case  this  is  interpreted  as  turbulence  doing  work  to 
mix  the  warmer  air  above  the  inversion  with  cooler  boundary-layer  air  below.  Thus,  the 
extra  —Bags  provided  by  the  standard  Smagorinsky  model  is  the  cause  of  the  enhanced 
entrainment  rate  seen  with  this  model.  Excessive  entrainment  of  warm,  dry  air  across  the 
inversion  then  leads  to  the  observed  reduction  of  cloud  fraction  and  liquid  water  path. 

The  same  effect  is  not  seen  with  the  dynamic  models  which  damp  the  model  coefficient 
so  that  it  becomes  zero  just  below  the  inversion.  Consequently  there  is  no  overlap  region 
and  no  region  of  significant  negative  SGS  buoyancy  production.  By  coincidence,  the 
same  effect  is  also  achieved  when  no  SGS  model  is  used,  which  explains  the  similarity 
between  the  results  obtained  with  the  dynamic  models  and  no  SGS  model  for  large  scale 
parameters  such  as  entrainment  rate,  liquid  water  path  and  cloud  fraction. 

Figure  4 shows  vertical  profiles  for  the  second  and  third  moments  of  the  resolved  ver- 
tical velocity  component  compared  with  the  observational  measurements.  The  second 
moment  is  well  predicted  in  the  dynamic  model  simulations,  with  numerical  results  in 
excellent  agreement  with  observations.  The  shape  of  the  profiles  obtained  with  the  dy- 
namic models  is  also  typical  of  a well-mixed  boundary  layer  as  is  expected  for  nocturnal 
stratocumulus.  In  contrast,  the  simulation  using  the  standard  Smagorinsky  model  does 
not  agree  as  well  with  the  observations  and  gives  a two-peaked  profile.  Such  a profile 
indicates  a decoupling  of  the  cloud  layer  and  sub-cloud  layer  which  is  not  apparent  in 
the  observational  results. 
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The  differences  in  profiles  of  the  third  moment  are  even  more  striking.  The  standard 
Smagorinsky  model  gives  reasonable  results  in  the  surface  layer,  but  does  not  give  the 
negative  values  observed  in  the  cloud  layer.  The  dynamic  models,  on  the  other  hand,  give 
good  agreement  with  the  observations  throughout  the  boundary  layer,  although,  close  to 
the  surface  they  appear  to  under-predict  the  magnitude  of  this  statistic  somewhat.  More 
work  is  required  on  the  surface  layer  model  to  improve  this. 


4.  Future  work 

The  results  presented  in  this  paper  clearly  demonstrate  the  utility  of  dynamic  turbu- 
lence models  for  atmospheric  simulations.  There  is,  however,  much  work  still  to  be  done. 
At  present  we  are  using  a second-order  accurate  scheme  for  the  spatial-discretisation. 
Consequently  the  numerical  errors  are  of  a similar  magnitude  to  the  terms  associated 
with  the  subgrid  scale  models.  The  authors  are  currently  working  on  a fourth-order 
scheme  to  resolve  this  issue.  More  work  is  also  required  in  refining  the  surface  layer 
model,  especially  where  it  is  applied  in  combination  with  the  dynamic  models.  As  dis- 
cussed above,  dynamic  models  tend  to  underpredict  the  dissipation  in  the  surface  layer. 
This  is  because  as  the  surface  is  approached  the  energetic  eddies  are  no  longer  resolved 
and  the  assumptions  underlying  the  dynamic  procedure  break  down.  While  the  current 
surface  layer  model  helps  to  overcome  these  issues,  the  surface  layer  velocity  and  scalar 
profiles  still  deviate  significantly  from  profiles  predicted  by  similarity  theory. 
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